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ABSTRACT 



This thesis analyzes the sway, yaw, and roll dynamic stability of neutrally 
buoyant submersibles. Utilizing the hydrodynamic coefficients for a Mark 
IX Swimmer Delivery Vehicle (SDV) as a base-line model, the linearized 
equations of motion for the decoupled steering and roll systems are 
compared to the coupled system. Two different configurations of 
hydrodynamic coefficients are considered along with the effects of varying 
the vertical (Zg) and longitudinal (Xg) centers of gravity of the vehicle while 
the longitudinal center of buoyancy (Xb) is held constant. Results 
demonstrate the significant effects on stability of coupling the steering and 
roll equations of motion, and the importance of Zg and Xg selection in 
minimizing those effects while retaining stability. Perturbation analysis 
results confirm the essential dependence of the linearized coupled 
equations on the separation of Xg and Xb- 
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I. D4TRODUCTION 



A. GENERAL 

As the missions for snbmersibles become more complex and 
demanding, the requirement for a highly stable platform becomes 
increasingly important so the operator(s) can concentrate on matters other 
than station-keeping- Submersible simulators have not been employed to 
any great extent thus far, since the actual system is relatively inexpensive 
and the dynamics are usually veiy non-linear and difficult to model. The 
analysis of a submersible can be significantly more complex than the 
analysis of a conventional submarine or aircraft, since the presence of 
ancillary equipment such as manipulators, video devices, and tethers 
introduce extra cross-coupling terms usually absent in other, more 
symmetric, vehicle dynamic analyses. [Ref. 1] Additionally, all 
mathematical models include simpiif^ng assumptions and errors in the 
model’s hydrodynamic coefficients. 

Submersibles typically have a variety of complex dynamic interactions 
that can severely inhibit maneuverability and control performance. The 
goal of this thesis is to present an understanding of the coupling effects on 
straight line motion stability in the horizontal plane using a linearized 
model, and the primary means of minimizing these effects. Development of 
the mathematical models for both the coupled and uncoupled maneuvering 
and roll equations of motion is presented in Chapter IT. Utilizing two 
different configurations of hydrodynamic coefficients, the degree of stability, 
regions of stability, and linear simulations for the coupled and uncoupled 
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systems are presented in Chapter III. A Hamming method nonlinear 
simulation [Ref. 2], which is similar to a fourth order Runge-Kutta 
integration technique, is conducted on both configurations to compare with 
the results obtained for the linearized models. Chapter IV develops a 
perturbation analysis to demonstrate the strong degree of dependence on 
the separation between the longitudinal centers of buoyancy and gravity to 
the solution of the linearized coupled equations of motion. Chapter V 
summarizes the results and provides recommendations for future 
submersible modelling research. Appendix A contains the computer 
programs utilized for the linear and nonlinear simulations. 

B. PARAMETER DEFINITION 

The values for the hydrodynamic derivatives and vehicle dimensions 
are from Smith, Crane, and Summey [Ref. 3], with the following exceptions: 

• Yr - the force in sway due to a change in yaw rate ^ 

• Nv - the moment due to a change in sway velocity. 

These two coefficients were modified to produce two different models that 
would have one eigenvalue change sign for a reasonable range of 
longitudinal and vertical centers of gravity. A comparison between the 
actual non-dimensional values and those used in this thesis is as follows: 

_Yr Nv_ 

Configuration ‘A’ -3.500E-02 -1.484E-03 

Configuration ‘B’ -5.940E-02 -1.484E-03 

Actual SDV 



2.970E-02 



-7.420E-03 



The effects of changing these coefficients are illustrated in Figures 1 and 2 
on the following pages. Additionally, the analysis presented herein is 
conducted in dimensional form; hence the nominal forward longitudinal 
velocity 'U’ appears in the equations of motion. All calculations and 
simulations utilize a value of 5 ft/sec for ‘U’. The coordinate system 
convention is the standard body-fixed, right -hand orthogonal axis system 
emplo}dng the Euler angle approach. 



1. Variables 



X, y, z 



u, V, w 



p» q.T 
X, Y, Z 
K, M,N 
V, 0, <t> 



Xg,Yg,Zg 

Xb. Yb, Zb 



Distances along the body fixed principal axes. 

Translational velocity components of model relative to 
fluid along body axes. 

Rotational velocity components of model along body axes. 

Hydrodynamic force components along body axes. 

Hydrodynamic moment components along body axes. 

Yaw, pitch, and roll angles; positive values following the 
right-hand rule. 

Center of gravity coordinates along principal axes. 
Center of buoyancy coordinates along principal axes. 



Ixx> fyyj Izz 
Xnose> Xtail 



h(x), b(x) 



Moments of inertia about the principal axes. 

Distances from body center along the longitudinal axis 
used in the crossflow force and moment integrals. 

Values are located within the nonlinear computer 
simulation program in Appendix A. 

Model width and height values used in the crossflow 
force and moment integrals. Values are located 
within the nonlinear computer simulation program 
in Appendix A. 
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DEGREE OF STABILITY 

Figure 1: Effect on Degree of Stability for Change in Yr and Nv - Configuration ‘A’ (Zg = 0.10 ft) 
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II. STABILITY OF MOTION 
A. EQUATIONS OF MOTION 

The horizontal plane, nonlinear equations of motion for a submersible 
as developed by Smith, Crane, and Summey [Ref. 1] are shown below in 
Equations (2.1). 

m[v+ur-wp + Xg(pq + r)- Yg(p^ + r^ ) + Zg (qr - p)] = 

[ Yp p + Yj. f + Ypq pq + Yqr qr] + [ YyU v + Yv wV w + Y§ j.u ^ 5 r] + 

[YyV + Ypup+ Yrur+ Yyqvq + YwpWp + Y^rwr] + (W - B)cos 0 sin (|) - 

J*[CDy h(x)(v + xr)^ + Cdz b (x )( w - xq )^1 ^ dx 

V... Ucf(x) 



Sway: 

(2.1a) 



Yaw: l^r + (lyy-Ixx)PQ ' Ixy(P^-Q^) - 

Ixz(qr-p) + m[Xg(v+ur-wp) - Yg(u-vr + wq)] = 

[NpP + Nrf + Npqpq + Nqrqr] + 

[NyV + Npup + Nrur + Nyqvq + Nwpwp +NwrWi] + 

[Nyuv + Nywvw + N5rU^6r] + (XgW-XbB)cos 0 sincb + (YgW-YbB)sin0 - 

J fCDyh(x)(v + xr)^ + Cdz b(x)(w -xq)^] (x) dx + u^Nprop 

Ucf(x) 
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Roll: 

(2.1c) 



IxxP+ Q^*(Izz ■ lyy) + Ixy(p^"" Q) " Iyz(Q )" 
Ixz(PQ +r) + m[Yg(w -uq + vp)-Zg(v + ur- wp)] = 
[Kpp + Kff + Kpqpq+Kqrqr] + 

[KyV +KpUp+Krur+ Kyqvq + KwpWp+KwrWr] + 

[Kvuv+ Kvwvw] + (YgW - YhB)cos 0 cos - 
(ZgW-ZbB)cos 6 sin ({) +u^Kpj-op 



B. SIMPLIFICATIONS 

In order to obtain the linearized equations of motion about a level flight 
path, the following simplifications were utilized: 



• The translational velocity (w) and acceleration (w ) in the z-direction 
are zero. 

• The rotational velocity (q) and acceleration (q) in the y-direction are 
zero. 

• The acceleration in the longitudinal direction (li ) is zero. 

• The cross-products of inertia are zero by virtue of a body-centered 
coordinate system. 

• The submersible is neutrally buoyant so B = W. 

• The longitudinal center of buoyancy (Xb) and the vertical center of 
buoyancy (Zb) are located at the origin of the body-fixed coordinate 
system. 

• The lateral center of buoyancy (Yb) and the lateral center of gravity 
( Yg ) are located at the origin of the body-fixed coordinate system. 

• Dynamic stability analysis is performed with all controls fixed; 
hence, all forces and moments due to control surfaces are zero. 

• The angle of pitch (0) is sufficiently small for sin(0) to equal zero. 
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• From Smith, Crane, and Summey [Ref. 1] the propeller coefficients 
Kprop and Nprop are zero. 

C. COUPLED STABILITY EQUATIONS 

When the simplifying assumptions from the preceeding section are 
applied, the resulting linearized equations are: 

(2.2a) Sway: m[v +ur + Xg(r)- Zg(p)] = Yf 

(2.2b) Yaw: Izzf + mXg(v+ur) = Nf 

(2.2c) Roll: IxxP - mZg(v+ur) = Kf 

where the force and moment representations are given by: 

(2.3a) Sway Force: Yf = YyV + YyV + Yfi* + Yrr + Ypp + Ypp 

(2.3b) Yaw Moment: Nf = NyV +NyV + Ni-r + Nrr + Npp + Npp + (XgW-XbB)cp 

(2.3c) Roll Moment: Kf = K^v +KvV + Kj.r + Krr + Kpp + Kpp+ (ZgW-ZbB)cp 

Equations 2.2 and 2.3 may be combined to form the state space 
representation: 

Ax = Bx (2.4a) 

where x = [p (|) v r]^ (2.4b) 
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A = 



(2.4c) 



Ixx - 


0 


-(Kv+MZg) 


-Kr 


0 


1 


0 


0 


-(Yp+MZg) 


0 


M- Yv 


MXg-Yr 


. -Np 


0 


M)fe.Nv 


Izz “ Nr . 



'KpU 


ZbB-ZgW 


KvU 


U(MZg+Kr)‘ 


1 


0 


0 


0 


YpU 


0 


YvU 


U(Yr-M) 


.NpU 


XgW-XbB 


NvU 


U(Nr-MXg). 



(2.4d) 



The stability of the coupled, linearized system depends on the location of the 
four eigenvalues of det(B - XA) = 0, which has a characteristic equation of 
the form 

A>.4 + B>c3 + C?i2 + D?i + E = 0, (2.4e) 

where the coefficients are complicated permutations of the elements in 
matrices A and B. The values for A, B, C, D, and E are given in Equations 
(2.4f - 2.4j) in terms of lower case letters that represent entries in matrices A 
and B; they are explicitly defined in Table 1. 



A = a-(j-u - 1-r) + d-(u-h + o-l) - f-(r-h + o-j) (2.4f) 

B = e-(u-h + 0-1) - d-(u-i - w-h + o-m - p-1) - a-(j*w + k-u - 1-x - r-m) 

- b-(j-u - 1-r) + f-(r-i - x-h + o-k - p-j) - g-(r-h + o-j) (2.4g) 

C = a-(k-w - m-x) + b-(j-w + k-u - 1-x - r-m) + g-(r-i - x-h + o-k - p-j) 

+ c-(j-u - 1-r) - d-(q-l - w-i + m-p) - e-(u-i - w-h 4- o-m - p-1) 

+ f-(q-j - x-i + p-k) (2.4h) 
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D = g*(q*j - x-i + p-k) - f-q-k + d-q-m - e-(q-l - w-i + m-p) - b (k-w - m-x) 

- c-(j-w + k-u - 1-x - r-m) (2.4i) 

E= c*(k*w - m*x) + e-q-m - g-q-k (2.4j) 



TABLE 1. COEFFICIENT DESCRIPTOR VALUES 



a = Lfx-Kp 


b = KpU 


c = ZgW-ZbB 


d = MZg + Kv 1 


e = KvU 


II 


g = MZgU 


h = MZg + Yp 


II 


j = M-Yv 


k = YvU 


1 = MXg-Yr 


m = U(Yr-M) 


0 = Np 


p = NpU 


q = XbB-XgW 


r = MXg-Nv 


U = Izz - Nr 


w = U(Nr-MXg) 


X = NvU 



D. UNCOUPLED STABILITY EQUATIONS 

Using matrices A and B from the preceeding section (Equations 2.4c and 
2.4d), uncoupling the steering and roll equations is straightforward: 

STEERING EQUATIONS 



■ M-Y» MXg-Yf' 


V 




"YvU 


U(Yr-M) ■ 


V 


MXg-N^ la-Nf. 


_f _ 




.NvU 


U(Nr-MXg). 


_r_ 



ROLL EQUATIONS 



'Ixx-Kp 0“ 


p 




■RpU 


ZbB-ZgW 


'p" 




0 1_ 






1 


0 




(2.5b) 
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The characteristic equations for the uncoupled steering and roll equations 
are much simpler than that for the coupled system, and are given as 
Equations (2.5c and 2.5d) in terms of the hydrodynamic parameters rather 
than the descriptors used in the previous section. 

The steering characteristic equation form is: + BlX + Cl = 0, where 

Al = (M-Yv)azz-N^) - (MXg-Yr)(MXg-Nv) 

Bl = -[azz-Nr)(YvU)+(M-Yv)(Nr-MXg)(U)+ (2.5c) 

(Yr -M)(Nv -MXg)(U)+(Yr -MXgXNvU)] 

Cl = (YvU^XNr-MXg) - (NvU^)(Y,-M) 

The roll characteristic equation form is: + Br?l + Cr = 0, where 

Ar = ixx-Kp 

Br = -KpU (2.5d) 

Cr = ZgW 
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III. RESULTS AND DISCUSSION 



A. DEGREE OF STABILITY 

The effects of changing the longitudinal and vertical centers of gravity 
(while keeping the center of buoyancy fixed at the vehicle center) on the 
degree of stability for configuration ‘A’ is illustrated in Figure 3. Degree of 
stability as utilized in this thesis is defined as the maximum real value of all 
characteristic roots, with negative values indicating a stable situation. The 
degree of stability for the uncoupled system is represented by the dashed line. 
It can be seen that the critical value of Xg for which motions become unstable 
is clearly a function of the metacentric height Zg, whereas the uncoupled 
system predicts a constant value of Xg. Figure 4 displays the variation of the 
imaginary part of the root value for configuration ‘A’. 

Figures 5 and 6 are analogous to Figures 3 and 4 for configuration ‘B’; 
they show the degree of stability for varying Xg and Zg values. For this 
configuration, the degree of stability has a stronger dependence on the 
location of Xg and the value of Zg. For almost all positive values of Xg, the 
complex conjugate roots are increasing in value and eventually becoming 
positive; this indicates an oscillatory response that diverges when the degree 
of stability is positive. 

B. REGIONS OF STABILITY 

Figure 7 shows the region of stability for configuration ‘A’, with the 
imcoupled system represented by a dashed line. The uncoupled system 
predicts stability for all values of Xg greater than 0.18 ft, while the coupled 
system predicts an additional region enclosed by the triangular area to the 
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Xg (ft) 

Figure 3: Effects of Xg and Zg on Degree of Stability for Configuration A 
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Figure 4: Effects of Xg and Zg on Imaginary Part of Root for Configuration A 
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Xg (ft) 

Figure 5; Effects of Xg and Zg on Degree of Stability for Configuration 
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Xg (ft) 

Figure 6: Effects of Xg and 5^ on Imaginary Part of Root for Configuration ‘] 
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Xg (ft) 

Figure 7: Calculated Zg Location vs Xg for Configuration 



left of the dashed line. Figure 8 is analogous to Figure 7 for configuration ‘B’. 
The large discrepancies between the predicted regions of stability in this case 
occur for Zg values less than 0.045 ft. For small values of Zg, there are 
corresponding small regions of Xg where stability is predicted by the coupled 
equations but not by the uncoupled equations. The root values in this region 
are complex conjugates with very small real parts. Figures 9(a) and 9(b) 
illustrate the effect of co-locating Xb and Xg. This results in a degree of 
stability for the coupled system that is nearly identical to the uncoupled 
system. 

C. INTERPRETATION OF RESULTS 

It may be shown by applying Routh’s criterion [Ref. 4: pp. 211-218] that 
for a fourth order equation of the form + C?i^ -i- D?i -i- E to have 

roots with all negative real parts the following must apply: 

i. ) BCD - AD^ - EB^ > 0 

ii. ) E > 0. 

If the quantity 'E’ is less than zero, the system will become unstable and the 
resulting motion will be a simple divergence. If, however, the value of the 
quantity BCD - AD^ - EB^ is negative, the resulting instability will result in 
an oscillatory motion due to the presence of complex conjugate roots with 
positive real parts. 

For the coupled system of equations, the condition E = 0 yields: 

Zg = Xg -[Kv(Yr-M)/(YvNr - Nv(Yr-M))]. 

while the uncoupled system of equations reduces to a constant term 
expression for Xg: 

Xg = [(MYv)/(YvNr - Nv(YrM))]. 
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Xg (ft) 

Figure 8: Calculated Zg Location vs Xg for Configuration ‘ 
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Xg (ft) 

Figure 9(a): Effects of Xg and Zg on Degree of Stability when Xb Equals Xg for Configuration 



0.05 







Ainiavis ao aaaoHa 



21 



Xg (ft) 

Figure 9(b): Effects of Xg and Zg on Degree of Stability when Xh Equals Xg for Configuration ‘ 



This explains the differences in the regions of stability illustrated in Figure 7, 
since the above equations show a linear relationship between Zg and Xg for 
the coupled equations and a constant value for Xg for the uncoupled 
equations. When the value for Xg coincides with the value for Xb, the 
constant term ‘E’ in the coupled equations is reduced to that of the constant 
associated with the uncoupled equations, and the resulting predicted degree 
of stability no longer depends on the value of Zg. Substituting the coefficient 
values for ‘E’ serves to clearly illustrate the reduction: 

E = (ZgW)[(YvU^)(NrMXg) - (NvU^XYr-M)] + 

(KvU^XYr-M)(XbB-XgW) - (MZgUXYvUXXbB-XgW). 

When Xb = Xg and B = W for the neutrally buoyant case, the second and third 
terms are reduced to zero. When ‘E’ is then set equal to zero (the condition 
for determining where the real roots change from negative to positive), the 
dependence on the value Zg is removed and the expression for Xg is the 
constant given for the uncoupled equations. This reduction is the explanation 
for the appearance of Figures 9(a) and 9(b). When the longitudinal centers of 
buoyancy and gravity coincide for a neutrally buoyant vehicle, the degree of 
stability for the coupled system of equations covering all metacentric height 
values is equivalent to the uncoupled system of equations. 

A simple reduction of the equation resulting from Routh’s criterion to 
determine when a pair of complex conjugate roots crosses the zero axis is not 
easily accomplished. Figure 10 is presented as confirmation that the 
locations of Xg for which BCD - AD^ - EB^ = 0 matches the locations given 
graphically in Figure 5. 
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Xg (ft) 

Figure 10: Normalized Routh Criterion Values vs Xg for Varying Zg Values - Configuration 



D. LINEAR SIMULATIONS 

Figures 11 through 14 present comparisons of the coupled and uncoupled 
system responses for configuration ‘A’, For the stable cases Xg = 0.40 fl and 
for the unstable cases Xg = -0.20 ft, while Zg is 0.20 ft for both. The unstable 
coupled case (Figure 11) illustrates a simple divergence for both angle of drift 
and roll angle. This is expected since the roots are not complex conjugates. 
The unstable, uncoupled case accurately predicts the divergence in angle of 
drift, but the roll response is predicted to be stable. This may be explained by 
examining the uncoupled equation of motion in roll: 

b' - b' (Kpu)/(ixx “Kp) + ({) (ZgW)y(ixx -Kp) = o. 

Substituting values for configuration ‘A’ yields: 

b" + (})' (1.475) + ({) (0.720) = 0. 

The solution of this equation results in a natural frequency of 0.85 rad/sec 
and a damping factor of 0.869. This is an underdamped case, since both roots 
have negative real parts and are complex conjugates. Substituting values for 
configuration ‘B’ yields: 

b" + b' (1.475) + (j) (0.144) = 0. 

The solution for this case results in a natural frequency of 0.380 rad/sec and a 
damping factor of 1.941, which represents an overdamped situation. This 
also demonstrates that a vehicle’s natural frequency in roll may be increased 
by increasing the metacentric height. The effects of the underdamping may 
be seen in Figures 12 through 14. 

Figures 15 through 20 provide a comparison between the coupled and 
uncoupled equations for configuration ‘B’. The effects of the overdamping are 
evident in figures 15 through 17. This set of figures demonstrate the 
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TIME (seconds) 

Figure 11: Roll and Angle of Drift vs Time for Configuration ‘A’ Coupled System 
ZfT = 0.20 ft Xg = -0.20 ft Initial Roll Angle = 1.0 degree 
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Figure 12: Roll and Angle of Drift vs Time for Configuration ‘A’ Coupled System 
Zg = 0.20 fl Xg = 0.40 ft Initial Roll Angle = 1.0 degree 
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TIME (seconds) 

Figure 13: Roll and Angle of Drift vs Time for Configuration ‘A’ Uncoupled System 

Z„ = 0.20 ft X„ = -0.20 ft 
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TIME (seconds) 

Figure 14: Roll and Angle of Drift vs Time for Configuration ‘A’ Uncoupled System 

Z„ = 0.20ft X„ = 0.40 ft 
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Figure 15: Roll and Angle of Drift vs Time for Configuration ‘B’ Uncoupled System 

Z„ = 0.04 a X„ = 0.20 ft 
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TIME (seconds) 

Figure 16: Roll and Angle of Drift vs Time for Configuration ‘B’ Uncoupled System 

Z„ = 0.04 a X. = LOO ft 
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Figure 17: Roll and Angle of Drift vs Time for Configuration ‘B’ Uncoupled System 

Zg = 0.04ft X„= 1.50 ft 
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Figure 18: Roll and Angle of Drift vs Time for Configuration ‘B’ Coupled System 
Zg = 0.04 fl Xg = 0.20 ft Initial Roll Angle = 1.0 degree 
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Figure 19: Roll and Angle of Drift vs Time for Configuration ‘B’ Coupled System 
Zp = 0.04 ft Xg = 1.00 fl Initial Roll Angle = 1.0 degree 
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TIME (seconds) 

Figure 20: Roll and Angle of Drift vs Time for Configuration ‘B’ Coupled System 
Zg = 0.04 ft Xg = 1.50 ft Initial Roll Angle = 1.0 degree 



magnitude of the discrepancies when the coupling effects are not considered. 
The comparison below summarizes the results of the figures representing 
simulations for configuration ‘B’, where ‘C’ stands for coupled and ‘UC’ for 
uncoupled. 



£ UC Q UC c nc 



Xg 


0.20 


0.20 


1.00 


1.00 


1.50 


1.50 


Roll Stable 


Y 


Y 


N 


Y 


N 


Y 


Drift Stable 


Y 


N 


N 


N 


N 


Y 



As was seen before, the effects of the coupling on the system results in 
complex conjugate eigenvalues with positive real parts. Therefore, the larger 
values of Xg result in increasingly divergent oscillations instead of the 
stability predicted by the uncoupled system. 

Figure 21 is a three-dimensional presentation of the roll amplitude vs 
time for configuration ‘B’ as Xg varies from 0.15 ft to 1.50 ft. This mesh 
capability of MATLAB allows a comparative view of several solutions, and 
the behavior of the roll response for the coupled system is easier to discern. 

E. NON-LEVEAR SIMULATIONS 

In order to provide a measure of the accuracy of the results obtained 
utilizing the linearized equations of motion, a simulation program for the 
non-linear equations of motion was developed using Hamming’s method 
[Ref. 2]. Hamming’s method utilizes a Milne predictor and incorporates a 
modifier step prior to the correction step. The primary advantage of using 
Hamming’s method is that only two derivative function evaluations are 
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f'igUre 21l Mesh Plot of Roll Angle vs Time for Configuration ‘B’ as Xg varies from 0.15 to 1.50 fl 



required per step rather than the four or more normally required by other 
popular methods. The local error is of the same order of magnitude (h^) as a 
more time consuming process such as a fourth order R\mge-Kutta, but the 
reduction in fimction evaluations results in a faster simulation. The formula 
is presented below, and may also be found in the nonlinear computer 
program simulation in Appendix A. 

HAMMING’S METHOD 
y(i+l)predicted = y(i-3) + (4h/3)[2«i) - fti-1) + 2fi;i-2)] 
y(i+l)modified — y(i+^)predicted * (112/121)[y(i)ppedicted ' 3^i)corrected] 
y(i+l)corrected = (l/8)[9y(i) - y(i-2) + 3h{fli+lXnodified + 2fi:i) - f(i-l))] 
y(i+l) = y(i+l)corrected + (9/121)[y(i+l)predicted ‘ 3^i+l)corrected] 

The first four values must be determined by another method; the Euler 
linear solution with a small step size proved sxifficient. 

Figure 22 represents the simulation for the unstable representation of 
configuration ‘A’. Rather than the exponentially increasing roll angle and the 
-90 angle of drift computed with the linear simulation, the nonlinear solution 
predicts an angle of drift that reaches -15 degrees, and then slowly diverges. 
The roll angle reaches a steady state value of approximately three degrees. 
Figures 23 and 24 are the nonlinear simulation results for stable 
configurations of ‘A’ and ‘B’, respectively. They are nearly identical to the 
results obtained using the linear simulation and displayed as Figures 12 and 
18. 

Figure 25 is the simulation for an unstable configuration ‘B’. Quite 
notable are the steady state roll angle and angle of drift after approximately 
250 seconds rather than the exponentially increasing divergence apparent in 
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TIME (seconds) 

Figure 22: Roll and Angle of Drift vs Time for Configuration ‘A’ - Nonlinear Simulation 
Zg = 0.20 fl Xg = -0.20 ft Initial Roll Angle = 1.0 degree 
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Figure 23: Roll and Angle of Drift vs Time for Configuration ‘A’ - Nonlinear Simulation 
Zg = 0.20 ft Xg = 0.40 ft Initial Roll Angle = 1.0 degree 
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TIME (seconds) 

Figure 24: Roll and Angle of Drift vs Time for Configuration ‘B’ - Nonlinear Simulation 
Zg = 0.04 ft Xg = 0.20 ft Initial Roll Angle = 1.0 degree 



angle Of drift 




c 

c 



o 

\Ti 

(N 



(S 



- m 



c 

o 



c 3 



c/; 

•o 

c 

o 

u 

o 

v: 

tl) 

2 

H 



__! O 

m 



:3 

g 

S 

cC 

O 

£ 

^ o 

g fe 

' 

p'i? 

C o> 

- 

o .£ 

V. .ti 

o 

H o 

03 O 
> 

II 

o>r 

Un 

o 

G) ^ 

'Hc'^ 

T3 II 
C bc 
03 t<j 

P 5 

u:: 

(N 

O 



(sasjSap) 310NV 



41 



c 




(S93JS3P) 310NV 



42 
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Figure 26: Roll Angle vs Time Comparison Between Linear and Nonlinear Simulations 

for Configuration ‘B’. Zg = 0.04 ft Xg = 1.00 ft Initial Roll Angle = 1.0 degree 
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Figure 27: Roll and Angle of Drift vs Time for Configuration ‘B’ - Enlarged View of Figure 25 

Z. = 0.04ft X.= 1.00 ft 
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Figure 28: Limit Cycle Phase Plane Trajectory for Roll Angle vs Angle of Drift - Configuration'' 

Zg = 0.04 ft X„ = 1.00 ft 



the corresponding linear simulation. Figure 26 provides a comparison 
between the linear and nonlinear solution for angle of roll. The enlarged 
view of the steady state region provided in Figure 27 better illustrates the 
constant limits of oscillation. By plotting the roll angle versus the angle of 
drift, an elliptical trajectory is apparent in Figure 28. This result is similar 
to the results obtained by Schmidt and Wright in their analysis of high 
performance aircraft wing-rock [Ref 5]. They postulate that a possible 
explanation for the limit cycle is the inertial coupling between a stable 
longitudinal and an unstable lateral mode. Similar results in the work are 
attributed to dynamic as well as hydrodynamic coupling. Figures 25 and 28 
indicate that the nonlinear interactions for the unstable conditions of 
configuration ‘B’ provide a significant amount of damping to the rolling 
motion. The limiting of the rolling motion accomplished by including the 
nonlinear terms then serves to limit the buildup of the angle of drift and the 
result is an eventual stable limit cycle. 

The ocean environment can be expected to introduce many combinations 
of disturbance forces, which may or may not be periodic. A preferred method 
for simulating many disturbances such as sensor noise, ocean current 
fluctuation, and vehicle acceleration fluctuations is to model them as ‘white 
noise’. Figure 29 is presented to demonstrate the effects of including 
constant, zero-mean disturbances in sway and yaw on the marginally stable 
system of configuration ‘B\ The disturbances are developed using MATLAB’s 
random number generator with a uniform di stribution. Che rnndom numbers 
are then scaled to simulate values that may be j^xpected. The variance of the 
sway and yaw acceleration disturbances, respectively, for Figure 29 are 0.003 
(fl/s^)^ and 0.002 (rad/s^)^. The resulting simulation bears little resemblance 
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to Figure 24, although the initial conditions and values for Xg and Zg are the 
same. With even these relatively minor disturbances acting on the system of 
configuration *B’, a rather large, non*symmetric oscillation in both roll and 
angle of drift is evident. Although the system is still stable, with the mean of 
both the roll angle and angle of drift equalling zero, a limit cycle similar to 
that of the imstable configuration (Figure 25) has developed. As the angle of 
drift fluctuates between positive three degrees and negative four degrees, the 
angle of roll varies between positive and negative two degrees. Increasing 
the scaling (which increases the variance) for the disturbances would be 
expected to increase the fluctuations imtil stability is lost. Similarly, it is 
true that the small disturbances acting in Figure 29 have a greatly reduced 
effect on the system of configuration ‘A’, which has greater stability. 
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Figure 29: Roll and Angle of Drift vs Time for Configuration ‘B’ with Constant Zero Mean Sway and 
Yaw Disturbances. Zg = 0.04 ft Xg = 0.20 ft Initial Roll Angle = 1.0 degree 



IV. PERTURBATION ANALYSIS 



A. BACKGROUND 

The previous section demonstrates that knowledge of the separation 
between longitudinal centers of buoyancy and gravity is critical in 
determining system stability. If this quantity is the dominant factor in a 
stability analysis, then an approximation for the degree of stability may be 
developed by applying perturbation theory. Although perturbation results 
are only approximations, their advantage over numerical methods is in 
illustrating the degree to which a solution depends on the variable(s) 
involved. From the fundamental theorem of perturbation theory [Ref 6], we 
seek a solution to the characteristic equation + CX^ + D>. + E = 0 of 

the form; 



n=oo 

X(e) = £ an e” 

n=0 

where e is the variable of interest and the solution approximates the 
numerical solution in the region where e is smailll. The constant coefficients 
(ao, ai, ..., an) axe all independent of £, and are alLtequal to zero for small e. 

B. PERTURBATION FORMULATION 

Recalling Equations 2.4f - 2.4j, which defimed the coefficients of the 
quartic characteristic equation, it would obviously be desirable to simplify the 
equations as much as possible prior to forimulating the perturbation 
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approximation. The variable of interest will be defined as: 

e = Xg - Xb (4.2) 

where the nominal value for Xb is zero and the perturbation will be 

performed around e = 0. By comparing Figures 30 and 31 it is clear that by 
allowing the hydrodjmamic coefficients Kv, Kr, Yp, Np, Yp, and to equal 

zero a very good approximate solution to Equation (2.4e) is obtained for both 
configurations ‘A’ and ‘B’. This simplification reduces the descriptor 
coefficients T, ‘i’, ‘o', and ‘p' from Table 1 to zero, which simplifies the 
coefficients of the coupled equations of motion to Equations (4.3a) - (4.3e) 
below. 

A = axx- Kp)[(M-Yv)(Izz-N,)-(MXg-Y,)(MXg-Nv)l + 

(MZg)2(I„-Nr) (4.3a) 

B = (KvU)(MZg)(I„-N,) + (MZg)2(Nr-MXg)(U)- 

(MZg)(MZg + Kr)(MXg - N J(U) - - Kp )[(M - YvXNr - MXg)(U) + 

(YvU)(I„ - Nr) - (NvU)(MXg - Yr) - (MXg -Nv)(Yr -M)(U)] - 
(KpU)[(M - Yv)(I^z - Nr ) - (MXg - Yr )(MXg - Nv )] 

C = (Ixx-Kp)[(YvU2)(Nr-MXg)-(NvU2)(Yr-M)] + 

(ZgWX(M-Yv)(I„-Nf)-(MXg -Yr)(MXg -Nv)] + 

(KpU)[(M - Yv)(Nr - MXg )(U) + (YvU)(I^r - Nr ) - 

(NvU)(MXg-Yr) - (MXg-Nv)(Yr-M)(U)] + 

(MZg)[(XgV,0(MXg -Y4-(NvU2)(MZg + Kr)+(KvU^)(Nr -MXg)] 



(4.3b) 



(4.3c) 



49 



D = -(ZgW)[(M-Yv)(Nr-MXg)(U) + (YvU)(I„-Ni.)- 
(N vUXMXg - Y^) - (MXg - N»)(Yr - M)(U)] - 

(KpU^)[(Yv)(N,-MXg)-(Nv)(Y,-M)] + (KvU)(XgW)(MXg-Y,) 
(MZg)(XgW)(U)(Y,-M) - (MZg+Kr)(U)(XgW)(M-Yv) 

E = (ZgW)[(YvU2)(Nr-MXg)-(NvU2)(Yr-M)] - (KvXgWU^XY^ -M) ■ 
(YvXgWU^XMZg+Kr) 

Recalling the definitions of the coefficients for the uncoupled systems: 

Al = (M-Yvlfizz-Ni.) - (MXg-YrXMXg-Nv) 

Bl = -[(lzz-Nr)(YvU) + (M-YvXNr-MXgXU) + 
dr - MXN V - MXgXU) + (Yr - MXg XNvU)] 

Cl = (YvU^XNr-MXg) - (NvU^XYr-M) 

Ar = Ixx-Kp BR = -KpU CR = ZgW 

Substituting in Equations (4.3a) - (4.3e) yields: 

A = ArAl + (MZg)^Izz-Nr) 

B = ArBl + BrAl + (MZg)2(N,-MXgXU) + 

(MZgXdzz -NpXKvU)-(MZg +KpXMXg -NvXU)] 

C = ArCl + CrAl + BrBl + (MZg)[(XgWXMXg-Yf) + 

(KvU^XNr - MXg) -(NvU^XMZg + K^)] 

D = CrBl + ClBr + (KvUXXgWXMXg-Y;.) - 

(MZg + K,)fXgW)(U)(M-Yv) - (MZgXXgWXU)(Yr-M) 

E = CrCl + (YvU‘‘)(XgWXMZg + Kr) - (KvXgWU^XYr -M) 



(4.3d) 

(4.3e) 



(4.4a) 

(4.4b) 

(4.4c) 

(4.4d) 

(4.4e) 
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Xg (ft) 

Figure 30: Comparison of Degree of Stability vs Xg with Six Cross Coupling 
Coefficients Neglected - Configuration ‘A’ 
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Figure 31: Comparison of Degree of Stability vs Xg with Six Cross Coupling 
Coefficients Neglected - Configuration ‘B’ 



In order to further simplify these coefficients, terms of order {Zgf and (Xg)^ 

will be neglected; the effects of doing this are small and may be seen in 
Figures 32 and 33. This reduces the coefficients to: 



where: 



A = 


ArAl 


(4.5a) 


B = 


ArBl + BrAl + a (Xg) + K1 


(4.5b) 


C = 


ArCl + CrAl + BrBl + P (Xg) + K2 


(4.5c) 


D = 


CrBl + ClBr + y (Xg) 


(4.5d) 


E = 


CrCl + 6 (Xg) 


(4.5e) 


a = 


■U%K,U 


(4.6a) 


P = 


-(MZg)[(WYf) + (MKvU^)] 


(4.6b) 


Y = 


(WUXMZg(Y«-Yr)+ KrYv-KrM-KvYf] 


(4.6c) 


5 = 


(WU^)[(Yv)(MZg+Kr) - (KvXYrM)] 


(4.6e) 


K1 


= (MZg)[(h^-N^)(KvU) + (NvKrU)] 


(4.6f) 


K2 = 


= (MZgU^)(NrKv - NvKr) 


(4.6g) 



Carrying through the computations results in the following expressions: 
(ArAl) + (ArBl+BrAl+K 1) ao^ + (ArCl+CrAl+BrBl+K 2) ao^ + 



(CrBl+ClBr) ao + CrCl = 0 



(4.7a) 



-(aap^ + Pap^ > yap -h 5) 

" (4A)ao®+3(B-oXg)ao^+2(C-pXg)ao + (D-7Xg) 



(4.7b) 



a2 



-ai^(C-PXg) - ai^ap(6A +B-gXg) - Klapai^ - 3aao^ai-2Paoai-Yai 
(4A)ao® + 3(B - aXg)ao^ + 2(C - pXg)ao + (D - yXg) 

(4.7c) 



53 



0.15 






o 








(S 




cn 




o 




o 


d 




d 




d 




d 




d 




d 


1 


d 


1 


d 



AiiTievisjoaH^oaa 



54 



Xg (ft) 

Figure 32: Comparison of Degree of Stability vs Xg Neglecting Higher Order Terms 
of Xg and Zg for Configuration ‘A’ 
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Xg (ft) 

Figure 33: Comparison of Degree of Stability vs Xg Neglecting Higher Order Terms 
of Xg and Zg for Configuration ‘B’ 



The characteristic equation + DX. + E = 0 may now be 

expressed in terms of a second order perturbation which depends on a single 
variable e as: 

Xfe) = ao + ai e + a 2 (4.8) 

Figures 34 through 37 depict the results for both first and second order 
perturbation analysis about Xg close to zero. Results for both configurations 

‘A’ and ‘B’ demonstrate that this method for representing the degree of 
stability in the vicinity of a nominal (Xg - Xb) is quite satisfactory, provided 

the roots have no complex conjugates. In the region where the roots become 
complex, the problem is no longer a regular perturbation, and methods to 
solve the singular perturbation must be developed. 
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Figure 34: First Order Perturbation Approximation for Configuration ‘A’ 
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Figure 35: First Order Perturbation Approximation for Configuration ‘1 
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Figure 36: Second Order Perturbation Approximation for Configuration ‘A’ 
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Figure 37: Second Order Perturbation Approximation for Configuration ‘1 



V. CONCLUSIONS AND RECOMMENDATIONS 



• When considering an inherently stable situation for an analysis of sway, 
yaw, and roll response, the linear simulations compare favorably with 
the nonlinear simulations. For a marginally stable or imstable system, 
only the nonlinear simulation c£m predict the existence and extent of any 
non-zero steady states or limit cycles that may occur. 

• Obviously, in the design of a submersible, careful attention must be paid 
to selection of both the metacentric height and the separation between the 
longitudinal centers of gravity and buoyancy. Instability associated with 
the dynamic coupling effects can be minimized by increasing the 
metacentric height. The uncoupled system stability predictions are not 
reliable when there is separation between Xg and Xb. 

• It can be concluded from the analysis that the coupling of roll into sway 
and yaw for the linearized equations of motion is not very significant 
when Xg equals Xb. 

• Recommendations for future modelling research include an expansion 
to incorporate coupling effects between the vertical and horizontal planes. 
It is well known that the coupling between translation and attitude is a 
severe limitation in the design of submersibles. 

• The effects of varying hydrodynamic derivatives and initial conditions, as 
well as incorporating disturbances in the nonlinear simulations deserve 
more attention. Additionally, the complexity of the nonlinear system 
resulting from including the vertical plane is greatly increased. More 
studies of the nonlinear system resulting from this union are also 
required. 
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APPENDIX A 



NON LINEAR SIMULATION PROGRAM: 



9): ^ 9|: 9|c 9): :|c 4: 9|c :|c 4: 4: :1c :)c D|c 9): 9|c 9|c ^ 9): * 9|( 9|c 9|c 9|c 9|c 9|c 9|c 9|c 9|( * * :|c :|c 9|c :|( 9): 9): 9): ^ 9|c 9): 9|c:tc 9|( 9|( 9|( 9|( 9|( * ^ 9|c « 3|( 9|; 9|( 9|( 9|( 9)c 9|c ;|c 9|c 9|c 9|( 9)( 

% THIS PROGRAM USES HAMMING’S METHOD TO SOLVE A SYSTEM 
% OF NONLINEAR EQUATIONS. SIMILAR TO A FOURTH ORDER 
% RUNGE-KUTTA IN TERMS OF LOCAL ERROR, IT REQUIRES ONLY 
% TWO FUNCTION CALCULATIONS PER STEP VICE FOUR. 

function[v,phi] = nonlinsim_hamming(UO,Xg,Zg) 

% VEHICLE PARAMETERS FOLLOW 



W = 12000,0; Ixx = 1760.0; lyy = 9450.0; Izz=10700.0; L = 17.425; 
RHO=1.94; G=32.2; M = W/G; BUOY = W; Zb=0.0; Xb = 0.0; 

ZZ = 0.5*RHO*L'^2; 

% SWAY HYDRODYNAMIC COEFFICIENTS 



Ypdot = 1.270E-04*ZZ*L'^2; 


Yrdot 


= 1.240E-03*ZZ*L^2; 


Ypq = 4.125E-03*ZZ*L^2; 


Yqr 


= -6.510E-03*ZZ*L'^2; 


Yvdot = -5.550E-02*ZZ*L; 


Yp 


= 3.055E-03*ZZ*L; 


Yvq = 2.360E-02*ZZ*L; 


Ywp 


= 2.350E-01*ZZ*L; 


Ywr = -1.880E-02*ZZ*L; 


Yv 


= -9.310E-02*ZZ; 


Yvw = 6.840E-02*ZZ; 







% NOMINAL VALUE FOR Yr = 2.970E-02*ZZ*L 
% CONFIGURATION ‘A’: Yr = -3.500E-02*ZZ*L 
% CONFIGURATION ‘B’: Yr = -5.940E-02*ZZ*L 



CC = input(‘enter model configuration choice; either 1 for A or 2 for B’) 
if CC<2 

Yr = -3.500e-02*ZZ*L; 
else 

Yr = -5.940e-02*ZZ*L; 
end 



% ROLL HYDRODYNAMIC COEFFICIENTS 



Kpdot = -1.01E-03*ZZ*L^3; 


Krdot 


= 


Kpq = -6.93E-05*ZZ*L''3; 


Kqr 


= 


Kvdot = 1.27E-04*ZZ*L''2; 


Kp 


= 


Kr = -8.41E-04*ZZ*L''2; 


Kvq 


= 


Kwp = -1.27E-04*ZZ*L''2; 


Kwr 


= 


Kv = 3.055E-03*ZZ*L; 


Kvw 


= 



^3.37E-05*ZZ*L^3; 

1.68E-02*ZZ*L^3; 

-1.10E-02*ZZ*L'^2; 

-5.115E-03*ZZ*L^2; 

1.39E-02*ZZ*L'^2; 

-1.87E-01*ZZ*L; 
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% YAW HYDRODYNAMIC COEFFICIENTS 



Npdot = -3.370E-05*ZZ*L^3; 
Npq = -2.110E-02*ZZ*L'"3; 
Nvdot = 1.240E-03*ZZ*L^2; 
Nr = -1.640E-02*ZZ*L^2; 
Nwp = -1.750E-02*ZZ*L''2; 
Nv = -1.484E-02*ZZ*L; 



Nrdot = -3.400E-03*ZZ*L^3; 

Nqr = 2.750E-03*ZZ*L^3; 

Np = -8.405E-04*ZZ*L'^2; 

Nvq = -9.990E-03*ZZ*L'^2; 

Nwr = 7.350E-03*ZZ*L'^2; 

Nvw = -2.670E-02*ZZ*L; 



% THE FOLLOWING ARE USED FOR EVALUATING THE INTEGRALS 
% IN THE SWAY AND YAW EQUATIONS OF MOTION. X IS THE 
% DISTANCE IN FEET ALONG THE LONGITUDINAL AXIS AND ‘HH’ 
% IS THE VEHICLE HEIGHT. ALL VALUES ARE IN FEET. 



X(l)=- 105.9/12; 

X(4)=-94.3/12; 

X(7)=-66.3/12; 

X(10)=79.2/12; 

X(13)=91.2/12; 

X(16)=101.2/12; 

HH(1)=0.0/12; 

HH(4)=13.96/12; 

HH(7)=29.36/12; 

HH(10)=30.00/12; 

HH(13)=21.44/12; 

HH(16)=9.12/12; 



X(2)=-104.3/12; 

X(5)=-87.3/12; 

X(8)=-55.8/12; 

X(ll)=83.2/12; 

X(14)=95.2/12; 

X(17)=102.1/12; 

HH(2)=2.28/12; 

HH(5)=19.76/12; 

HH(8)=31.85/12; 

HH(11)=27.84/12; 

HH(14)=17.12/12; 

HH(17)=6.72/12; 



X(3)=-99.3/12; 

X(6)=-76.8/12; 

X(9)=72.7/12; 

X(12)=87.2/12; 

X(15)=99.2/12; 

X(18)=103.2/12; 

HH(3)=8.24/12; 

HH(6)=25.1/12; 

HH(9)=31.85/12; 

HH(12)=25.12/12; 

HH(15)=12.0/12; 

HH(18)=0.00/12; 



% THE INITIAL CONDITIONS ARE SET PRIOR TO INTEGRATION 
% BY CALLING A LINEAR INTEGRATION PROGRAM NAMED 
% ‘EULER’. THE ERRORS ASSOCIATED WITH USING LINEAR 
% SOLUTIONS FOR THE FIRST FIVE TIME INTERVAL STEPS ARE 
% SMALL AND DO NOT AFFECT THE NONLINEAR SOLUTIONS. 



[p,pdot,v,vdot,r,rdot,phi,phidot] = euler(dt); 

p(l:5)= p; 
pdot(l:5)= pdot; 
pdotmod(5)=pdot(5); 

v(l:5)= v; 
vdot(l:5)= vdot; 
vdotmod(5)=vdot(5); 
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^ ^ ^ ^ ^ ^ ^ 



r(l:5)= r; 
rdot(l:5)= rdot; 
rdotmod(5)=rdot(5); 



phi(l:5)= phi; 
phidot(l:5)= piildot; 
phidotmod(5)=phidot(5); 

% THE TIME INTERVAL IS ‘dt’, THE FINAL TIME IS ’tfinal’, AND A 
% VALUE CALLED ‘stopnumber’ IS SET TO ALLOW ACCESS TO THE 
% DATA AFTER THIS NUMBER OF ITERATIONS. THE VALUE FOR 
% ‘stopnumber’ MAY BE CHANGED WHILE THE KEYBOARD IS 
% ACTIVE TO ALLOW SUBSEQUENT PROGRAM INTERACTION. 

% RETURN TO THE PROGRAM IS ACHIEVED BY TYPING ‘return’ 
% FOLLOWED BY THE ‘enter’ KEY. 



dt = input(‘enter the time interval step size’) 

tfinal = input(‘enter the final time’) 

stopnumber = input(‘enter the value for stopnumber’) 

% THIS SECTION ALLOWS FOR RANDOM DISTURBANCES IN SWAY 
% AND YAW 

rand( uniform*) 

vdotdist = rand(l,(tfinal/dt)+l); 
vdotdist = 0.05* (vdotdist - mean(vdotdist)); 
rdotdist = rand(l,(tfinal/dt)+l); 
rdotdist = 0.04*(rdotdist - mean(rdotdist)); 

THIS SECTION PROVIDES FOR CONTINUATION OF SOLUTIONS. 
IF THE MATRICES BECOME VERY LARGE, IT IS RECOMMENDED 
THAT THE CURRENT VALUES BE SAVED AND A NEW 
SIMULATION BEGUN AS MATLAB SLOWS NOTICEABLY WITH 
LARGE MATRICES. 



vcorr(5)= 
pcorr(5)= 

% vpred(5)= 
% ppred(5)= 



rcorr(5)= 

phicorr(5)= 

rpred(5)= 

phipred(5)= 



% pdotmod(5)= 



rdotmod(5)= 
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% COMMENCE INTEGRATION 



for j = 6:(tfinal/dt) + 1; 

% THIS SECTION USES A HAMMING PREDICTOR-CORRECTOR TO 
OBTAIN VALUES 

vpred(j)=v(j-4)+(4*dt/3)*(2*vdot(j-l)-vdot(j-2)-h2*vdot(j-3)); 

rpred(j)=r(j-4)+(4*dt/3)*(2*rdot(j-l)-rdot(j-2)+2*rdot(j-3)); 

ppred(j)=p(j-4)+(4*dt/3)*(2*pdot(j-l)-pdot(j-2)+2*pdot(j-3)); 

phipred(j)=phi(j-4)+(4*dt/3)*(2*p(j-l)-p(j-2)+2*p(j-3)); 

ifj<7 



vmod(j) = 
rmod(j) = 
pmod(j) = 
phimod(j)= 
else 

vmod(j)= 

rmod(j)= 

pmod(j)= 

phimod(j)= 


vpred(j); 

rpred(j); 

ppredCj); 

phipred(j); 

vpred(j)-(112/121)*(vpred(j-l)-vcorr(j-l)); 

rpred(j)-(112/121)*(rpred(j-l)-rcorr(j-l)); 

ppred(j)-(112/121)*(ppred(j-l)-pcorr(j-l)); 

phipred(j)-(112/121)*(phipredy-l)-phicorr(j-l)); 



end 

tmplmod=vmod(j); tmp2mod=rmod(j); 

[SWAYMOD,YAWMOD]=crossflowintmod(tmplmod,tmp2mod,X,HH); 



vdotmod(j)= 


(l/(M-Yvdot))*((Ypdot+M*Zg)*pdotmod(j-l)+... 

Yp*UO*pmod(j)+(Yrdot-M*Xg)*rdotmod(j-l)-i-... 

Yv*UO*vmod(j)+(Yr-M)*UO*rmod(j)+SWAYMOD); 


vcorr(j)= 


0.125*(9*v(j-l)-v(j-3)+3*dt*(vdotmod(j)+2*vdot(j-l)-... 

vdot(j-2))); 


rdotmod(j)= 


(l/(Izz-Nrdot))*(Npdot*pdotmod(j-l)+..- 

(Nvdot-M*Xg)*vdotmod(j)-i-Np*UO*pmod(j)+... 

(Nr-M*Xg)*UO*rmod(j)+Nv*UO*vmod(j)+... 

Xg*W*sin(phimod(j))+YAWMOD); 


rcorr(j)= 


0.125*(9*r(j-l)-r(j-3)+3*dt*(rdotmod(j)+2*rdot(j-l)-... 

rdot(j-2))); 


pdotmod(j)= 


(l/(Ixx-Kpdot))*((Kvdot+M*Zg)*vdotmod(j)+... 

Krdot*rdotmod(j)+Kp*UO*pmodCL+... 

(Kr+M*Zg)*UO*rmod(j)+Kv*UO*vmod(j)-... 

Zg*W*sin(phimod(j))); 
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pcorr(j)= 0.125*(9*p(j-l)-p(j-3)+3*dt*(pdotmod(j)+2*pdot(j-l)-... 

pdot(j-2))); 

phicorr(j)= 0.125*(9*phi(j-l)-phi(j-3)+3*dt*(pmod(j)+2*p(j-l)-... 

p(j-2))); 

v(j)= vcorr(j )+(9/12 1 )*(vpred(j )-vcorr(j )); 
r(j)= rcorr(j)+(9/121)*(rpred(j)-rcorr(j)); 

p(j)= pcorr(j)+(9/121)*(ppred(j)-pcorr(j)); 

piii(j)= phicorr(j)+(9/121)*(phipred(j)-phicorr(j)); 

beta(j) = -atan(v(j)/5)*180/pi; 

% THIS SECTION HAS THE NON-LINEAR EQUATIONS OF MOTION. 
% REMEMBER TO REMOVE THE DISTURBANCES FROM THE SWAY 
% AND YAW DERIVATIVES IF THEY ARE NOT DESIRED. 

tmpl=v(j); tmp2=r(j); 

[SWAY, YAW] = crossflowint(tmpl,tmp2,X,HH); 

vdot(j)= (d/(M-Yvdot))*((Ypdot+M*Zg)*pdota-l)+Yp*UO*p(j) +... 
(Yrdot-M*Xg)*rdot(j-l)+Yv*UO*v(j)+(Yr-M)*UO*r(j)+... 
SWAY))-i-vdotdist(j); 

rdot(j)= ((l/(Izz-Nrdot))*(Npdot*pdot(j-l)-i-(Nvdot-M*Xg)*vdot(j)+-.. 
Np*UO*p(jH(Nr-M*Xg)*UO*r(j)+Nv*UO*v(j)-H... 
Xg*W*sin(phi(j))+YAW))+rdotdist(j); 

pdot(j)= (l/(Ixx-Kpdot))*((Kvdot-i-M*Zg)*vdot(j)-i-Krdot*rdot(j)-i-... 
Kp*UO*p(j)-i-(Kr-i-M*Zg)*UO*r(j)+Kv*UO*v(j)-... 
Zg*W*sin(phi(j))); 

ifj == stopnumber 
keyboard 
end 

end 

return 
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LINEAR SIMULATION PROGRAM: 



^^^♦♦^♦♦♦♦’f:****’*:’*:********************* ************************ ********** 

% THIS PROGRAM IS THE LINEAR EQUATION OF MOTION 
% INTEGRATION EMPLOYING A SIMPLE EULER METHOD TO 
% OBTAIN STARTING VALUES FOR THE NONLINEAR SIMULATION. 

function[p,pdot,v,vdot,r,rdot,phi,phidot] = euler(dt) 

% THE MATRIX ‘A’ CORRESPONDS TO THE DAMPING MATRIX AND 
% MATRIX ‘B’ CORRESPONDS TO THE MASS MATRIX. 



A(l,2) = -((Zg*WHZb*B)); 
A(l,4) = (Kr+(M*Zg))*U0; 
A(2,2) = 0.0; 

A(2,4) = 0.0; 

A(3,2) = 0.0; 

A(3,4) = (Yr-M)*U0; 
A(4,2) = ((Xg*W)-(Xb*B)); 
A(4,4) = (Nr-(M*Xg))*U0; 



A(l,l) = Kp*U0; 

A(l,3) = Kv*U0; 
A(2,l)=1.0; 

A(2,3) = 0.0; 

A(3,l) = Yp*U0; 

A(3,3) = Yv*U0; 

A(4,l) = Np*U0; 

A(4,3) = Nv*U0; 

B(l,l) = Ixx-Kpdot; 

B(l,3) = -(Kvdot+(M*Zg)); 
B(2,l) = 0.0; 

B(2,3) = 0.0; 

B(3,l) = -(Ypdot+(M*Zg)); 
B(3,3) = M-Yvdot; 

B(4,l) = -Npdot; 

B(4,3) = (M*Xg)-Nvdot; 

C = [inv(B)*A]; 



B(l,2) = 0.0; 

BQ,4) = -Krdot; 

B(2,2) = 1.0; 

B(2,4) = 0.0; 

B(3,2) = 0.0; 

B(3,4) = (M*Xg)-Yrdot; 
B(4,2) = 0.0; 

B(4,4) = Izz-Nrdot; 



% INITIAL CONDITIONS 



timed) = 0.0; 
p(l) = 0.0; 
phi(l) = 1.0*pi/180; 
psi(l) = 0.0; 
v(l) = 0.0; 

Til) = 0.0; 
xd) = 0.0; 
xdotd) = 0.0; 
tfinal = 1; 



pdotd) = 0.0; 
phidot(l) = 0.0; 
psidot(l) = 0.0; 
vdotd) = 0.0; 
rdotd) = 0.0; 
y(l) = 0.0; 
ydotd) = 0.0; 
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% THIS SECTION ALLOWS FOR RANDOM DISTURBANCES IN SWAY 
% AND YAW FOR THE LINEAR EQUATIONS AS WELL. MATLAB 
% ALLOWS FOR EITHER A UNIFORM OR NORMAL DISTRIBUTION 
% OF RANDOM NUMBER. 

rand(’imiforin') 

vdotdist = rand(l,(tfinal/deltat)+l); 
vdotdist = 0.05*(vdotdist - mean(vdotdist)); 
rdotdist = rand(l,(tfinal/deltat)+l); 
rdotdist = 0.04*(rdotdist - mean(rdotdist)); 

% COMMENCE EULER INTEGRATION 

for j = 2:(tfinal/deltat) + 1; 

P(j) = pG“ 1) + pdot(j-l)*deItat; 

phi(j) = phi(j-l) + phidoty-l)*deItat; 
psi(j) = psi(j-l) + psidot(j-l)*deItat; 
vG) = vG-1) + vdotG-l)*deItat; 

betaG) = -atan(vG)AJO); 
rG) = rG-1) + rdotG-l)*deItat; 

xG) = xG-1) + xdotG'l)*deItat; 

yG) = yG-1) + ydotG-l)*deltat; 



pdotG) = 

phidotG)= 

vdotG) = 

rdotG) = 

psidotG) = 

xdotG) = 

ydotQ) = 

timeG) = 

end 

return 



C(1,1)*pG) + C(l,2)*phiG) + C(l,3)*vG) + C(l,4)Mj) 
C(2,1 )*pG) + C(2,2)*phiG) + C(2,3)*vG) + C(2,4)*rG) 
C(3,1 )*pG) + C(3,2)*phiG) + C(3,3)*vG) + C(3,4)*rG) 
C(4,1 )*pG) + C(4,2)*phiG) + C(4,3)*vG) + C(4,4)*rG) 
rG)*cos(phiG)); 

U0*cos(psiG)) - vG)*sin(psiG))*cos(phiG)); 
UO*sin(psiG)) + vG)*cos(psiG))*cos(phiG)); 
G*deltat) - del tat; 
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CROSSFLOW INTEGRAL PROGRAM: 



♦ s|e *4: :i|c * :<c ♦ :iie :iic :iic :iic :<e sje sie sie :jc :jc * s|e :je + :<c sie sK * :je * :iic sic sie sie :iic :<c sie sie sie :iic s|e :je s|c s|e s|c sje :<c :<e :iic :iic sie sie sic * :jc 4: sit ^ sjc sic sj5 ;jc * 

% THIS IS A NUMERICAL INTEGRATION PROGRAM TO CALCULATE 
% THE DRAG FORCES IN THE HORIZONTAL PLANE UTILIZING THE 
% TRAPEZOIDAL RULE. THE CALL TO ‘crossflowintmod’ IS 
% IDENTICAL, ONLY USING DIFFERENT VARIABLE NAMES AND 
% VALUES FOR THE MODIFICATION PORTION OF HAMMING’S 
% METHOD. 



function! SWAY, YAW] =crossflowint(tmpl,tmp2,X,HH) 

RHO = 1.94; CDy = 0.35; SWAY = 0.0; YAW = 0.0; 

V = tmpl; r = tmp2; 



fork =1:18; 
if abs(v+X(k)*r) < le-6 
UCF(k) = 0.0; 
clso 

UCF(k) = (v+X(k)*r)/(abs(v+X(k)*r)); 
end 

CFLOW(k) = CDy*HH(k)*((v+X(k)*r)^2); 

SWAYl(k) = CFLOW(k)*UCF(k); 

YAWl(k) = CFLOW(k)*X(k)*UCF(k); 
end 

ii = 1:17; 

SWAY = -0.25*RHO*sum((SWAYl(jj)+SWAYl(jj+l)).*(X(jj+l)-X(jj))); 
YAW = -0.25*RHO*sum((YAWl(jj)+YAWl(jj+l)).*(X(jj+l)-X(jj))); 
return 



EIGENVALUE CALCULATION PROGRAM: 



% THIS PROGRAM DETERMINES THE FOUR ROOTS FOR THE ROLL 
% COUPLED EQUATIONS OF MOTION OVER A RANGE OF Xg. 

function[all_roots,xg] = findroots(UO,Xgmin,stepXg,Xgmax) 

W = 12000.0; Ixx = 1760.0; lyy = 9450.0; Izz=10700.0; 

L= 17.425; RHO=1.94; G=32.2; M = W/G; 

BUOY = W; Zb=0.0; 

ZZ = 0.5*RHO*L'^2; 
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% SWAY HYDRODYNAMIC COEFFICIENTS 

Ypdot= 1.27E-04*ZZ*L''2; Yrdot= 1.24E-03*ZZ*L''2; 

Yvdot = -5.55E-02*ZZ*L; Yp = 3.055E-03*ZZ*L; 

Yv = -9.31E-02*ZZ; 

CC = inputCenter model configuration choice; either 1 for A or 2 for B’) 
if CC<2 

Yr = -3.500e-02*ZZ*L; 
g1s0 

Yr = -5.940e-02*ZZ*L; 
end 



% ROLL mDRODYNAMIC COEFFICIENTS 

Kpdot = -1.01E-03*ZZ*L^3; Krdot = -3.37E-05*ZZ*L^3; 

Kvdot = 1.27E-04*ZZ*L^2; = -1.10E-02*ZZ*L'^2; 

Kr = -8.41E-04*ZZ*L^2; Kv = 3.055E-03*ZZ*L; 

% YAW HYDRODYNAMIC COEFFICIENTS 

Npdot = -3.370E-05*ZZ*L^3; Nrdot = -3.400E-03*ZZ*L^3; 
Nvdot = 1.240E-03*ZZ*L^2; Np = -8.405E-04*ZZ*L'^2; 
Nr = -1.640E-02*ZZ*L^2; Nv = -1.484E-02*ZZ*L; 



rows=(abs(Xgmax-Xgmin)/stepXg)+ 1; 
all_roots = zeros(rows,4); xg = zeros(rows,l); 



Zg = inputCenter value for Zg ) 



a = (Ixx-Kpdot); 
d = (M*Zg)+Kvdot; 
g = (M*Zg*U0) + (Kr*U0); 
i = Yp*U0; 

0 = Npdot; 

m = (Yr*U0) - (M*U0); 



b = (Kp*U0); c = (Zg*W)-(Zb*BUOY) 

e = Kv*U0; f= Krdot; 

h = Ypdot + (M*Zg); 
j = M - Yvdot; k = Yv*U0; 
u = Izz - Nrdot; p = Np*U0; 
x = Nv*U0; 



Xg = Xgmin:stepXg:Xgmax; 

Xb = input(‘enter either 0.0 or Xg for value of Xb); 

1 = (M*Xg) - Yrdot; q = (Xb.*BUOY)-(Xg*W); 

r = (M*Xg)-Nvdot; w = (Nr*U0)-(M*Xg*U0); 



A = (((a).*(j'*=u-l.*r)+(d).*(u*h+o.*l)-(f).*(r.’''h-rO*j))); 

B = ((e).*(u*h+o.*l)-(d).*(u*i-w.*h+o.*m-p.*l)-(a).*... 
(j.*w+k*u-l.*x-r.*m)-(b).*(j*u-l.*r)+(f).*(r.'^i-x... 
.*h+o*k-p*j)-(g).*(r.*h+o*j)); 
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C = ((a).*(k.*w-m.*x)+(b).*(j.*w+k*u-L*x-r.*m)+(c).*... 
o.*m-p.*l)+(f).*(q.*j-x*i+p*k)+(g).*(r.*i-x-*h+o*k... 



D = ((g).*(q.* 3 -x.*i+p*kXf-*q*k)+(d.*q.*m>-(e)-*(q.*l-. 

-w.*i+m.*pHc).*(j.*w+k*u-l.*x-r.*in)-(b).*(k.*w-in.*x)); 

E = ((c).*Ck.*w-m.*x)+(e.*q.*in)-(g.*q.*k)); 

z = [A' B’ C’ D’ E’]; 

j = Irrows; poly(3,.*)=z(j,*.); rts = roots(poly(ij,:))’; 

for jj = l:rows; 
all_roots(jj ,1:4) = eval(rts)’; 
end 

xg=Xg’; 

end; 

return 



FIRST AND SECOND ORDER PERTURBATION 
APPROXIMATION COMPUTER PROGRAM: 



% THIS PROGRAM DEVELOPS THE 1ST & 2ND ORDER PERTURBATION 
% APPROXIMATION SOLUTIONS FROM SECTION IV 

function[xlocation,perturbrootl,perturbroot2,xg] =perturbanalysis(UO) 

W = 12000.0; Ixx = 1760.0; lyy = 9450.0; Izz=10700.0; 

L= 17.425; RHO=1.94; G=32.2; M = W/G; 

BUOY = W; Zb=0.0; Xb = 0.0; 

ZZ = 0.5*EHO*L^2; 

% SWAY HYDRODYNAMIC COEFFICIENTS 

Ypdot= 0.0; Yrdot^ 1.24E-03*ZZ*L'^2; 

Yvdot = -5.55E-02*ZZ*L; Yp = 0.0; 

Yv =-9.31E-02*ZZ; 

CC = input(‘enter model configuration choice; either 1 for A or 2 for B’) 
if CC<2 

Yr = -3.500e-02^23*L; 
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6ls0 

Yr = -5.940e-02*ZZ*L; 
end 

% ROLL HYDRODYNAMIC COEFFICIENTS 
Kpdot = -1.01E-03*ZZ*L^3; Krdot - 0.0; 

Kvdot = 0.0; Kp = -1.10E-02*ZZ*L'^2; 

Kr = -8.41E-04*ZZ*L^2; Kv = 3.055E-03*ZZ*L; 

% YAW HYDRODYNAMIC COEFFICIENTS 

Npdot = 0.0; Nrdot = -3.400E-03^ZZ*L^3; 

Nvdot = 1.240E-03*ZZ*L^2; Np = 0.0; 

Nr = -1.640E-02*ZZ*L^2; Nv = -7,42E-03*EZ*L; 

Zg = inputCenter value for Zg') 

% THIS PERTURBATION PROGRAM SOLVES FOR THE SOLUTION 
% ABOUT Xg=0 

Xg = 0.0; 

index = 1; 



for Xg = -1.0:0.05:0.2; 

a = (Ixx-Kpdot); 
d = (M*Zg); 

g = (M*Zg*U0)-i-(Kr*U0); 
j = M - Yvdot; 
m= (Yr*U0) - (M*U0); 
q = (Xb*BUOY)-(Xg*W); 
w = (Nr*U0)-(M*Xg*U0); 



b = (I^*U0); 
e = Kv*U0; 
h = Ypdot ^ <M*Zg); 
k = Yv*U0; 

0 = Npdot; 
r = (M*Xg)-Nvdot; 

X = Nv^UO; 



c = (Zg*WHZb*BUOY); 
f= Krdot; 
i = Yp*U0; 
l = (M*Xg)-Yrdot; 
p = Np*U0; 
u = Izz - Nrdot; 



Cl = (k*w)-(x*xQ); 



A1 = (j*u)-(l*r); B1 = -((u*k)-»-(j*w)-(in*r)-(l*x)); 



Ar = a; Br = -b; Cr = c; 

K1 = d*u*e + d*Nvdot*Kr*U0; 

K2 = d*((Nr*Kv*U0^2MNv*Kr*U0^2)); 



alpha = -M*M*Zg*Kr*U0; 

beta = -d*((W*YrdotH(M*Kv*U0*U0)); 

gamma = (d*W*U0*(Yvdot-Yr))+(W*U0’^(Kr*Yvdot-Kr*M-Kv*Yrdot)); 
delta = g*Yv*W*U0 - m^V*W*U0; 
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% P Q R S &T SOLVE THE 4TH ORDER POLYNOMIAL FOR aO 



P = Ar*Al; 

Q = Ar*Bl + Br*Al +K1; 

R = Ar*Cl +Cr*Al +Br*Bl +K2; 

S = Cr*Bl + Cl*Br; 

T = Cr*Cl; 

A = p* 

B = Q+alpha*Xg; 

C = R+beta*Xg; 

D = S+gamma*Xg; 

E = T+delta*Xg; 

z = [A B C D E]; templ(index,l:4) = roots(z)’; 

xg(index,l) = Xg; index = index + 1; end 

zz = [P Q R S T]; temp2(l,l:4) = roots(zz)'; 

aO = temp2(l,4) 

% THIS PART SOLVES FOR al IN THE 1ST ORDER PERTURBATION 

numl = -(alpha*(a0^3)+beta*(a0'^2)+gamma*a0+delta); 
deni =(4^P*(a0'^3)+3*Q*(a0^2)+2*R*(a0)+S); 
al = numl/denl 

% THIS PART SOLVES FOR a2 IN THE 2ND ORDER PERTURBATION 

num2a= -(al^2*(6*P*aO*aO+Q*aO+R)+al*(gamma+2*beta*aO+... 

3*alpha*a0*a0)); 
a2 = num2a/denl 

index = 1; 

for X = -0.8:0.04:0.4; 



% PERTURBROOTl IS 1ST ORDER PERTURBATION APPROXIMATION 
perturbrootl(index) = aO + al*X; 

% PERTURBROOT2 IS 2ND ORDER PERTURBATION APPROXIMATION 
perturbroot2(index) = aO -i- al*X + a2*X*X; 

xlocation(index) = X; 
index = index+1; 
end 
return 
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